 module LowVOC

	using PlutoUI
	using ModelingToolkit
	using Unitful
	using Plots
	using Statistics
	using CSV
	using DataFrames
	using Markdown 

# Absorptive partitioning.

md"""
Kₒₘ is an equilibrium partition coefficient, which describes the partition between gas and aerosol phase for the compounds which can occurr in the form of particles and gases.
Naphtalene falls in the group of intermidiate volatile compounds in most of the cases which are more volatile than semivolatile intermidiate compounds. The partitioning between gas and aerosol phase can be calculated assuming pseudo-ideal solution and absorptive partitioning theory, which represents with the following equation:

Here we calculate Kₒₘ for semivolatile and intermediate volatile organic compounds, which represent by the following compounds: naphtalene, dibutyl phtalate, nanodecane. It corresponds to the equation #1 in the Pye et al. (2010).
"""

#naphtalene partial partition coefficient

#naphtalene partition coefficient (IVOC)


	R = 8.314 #unit = u"kg*m^2/(s^2*K*mol)" # gas constant
			#T is the temperature, unit = K
	Mᵢ = [128.1705; 78.11] #unit = ug/mol" #Molecular weight
	γᵢ = [2.3300; 4.54] #[unit = u"N*s^2/m^2"] # activity coefficient of compound i in aerosol phase
	T = [263.61, 267.98, 273.16, 278.22, 283.14, 288.01, 293.24, 293.25, 298.26, 303.29, 308.17, 313.24, 318.21, 323.14, 328.24, 333.39, 338.10, 343.06] #temperature range in K
	p = [0.23, 0.40, 0.74, 1.38, 2.41, 4.13, 6.93, 6.95, 11.35, 18.45, 28.95, 44.73, 68.82, 104.14, 158.41, 237.54, 340.76, 488.58] #vapor pressure in Pa
	Tp = hcat(T, p)
	C_sat = [1646.0, 20.0, 16.46, 0.2, 100000, 1.69, 270.0, 0.0001] #units saturation concentration
	species = ["SVOC1"; "SVOC2"; "O-SVOC1"; "O-SVOC2"; "IVOC"; "O-IVOC N1"; "O-IVOC N2"; "O-IVOC H1"]
	
	Kₒₘ = zeros(18)
	for n in 1:size(Tp,1), j in 1:size(Tp,2)
		Kₒₘ[n] = (R * Tp[n,1]) / (Mᵢ[1] * γᵢ[1] * Tp[j,2]*1000)
	end
	Kₒₘ #m3/ug
	
md """
Kₒₘ in m^3/ug is proportional to the gas constant, temperature, and reversely proportional to the molar mass, activity coefficient of the compound in aerosol phase, and partial vapor pressure. Here we calculate it using for loop for different temperatures. We can also plot it to see the relationship with the change of the temperature since the vapor pressure also changes with the temperature, the relationship is not that apparent. Now we can plot the equilibrium partition coefficient with the change of the temperature.
"""

plot(Kₒₘ, T, c=:blue, lab = :none, xlabel="Kₒₘ", ylabel="Temperature")

md """
Since the partitioning in under the highly dilute conditions found in the atmosphere, is strongly influenced by the ambient temperature. The plot gives the dependance of equilibrium air partitioning coefficient with the temperature.

SVOCs from all sources are assumed to be emitted as two semivolatile surrogate species, SVOC1 and SVOC2, in roughly equal fractions of 0.49 and 0.51. Under most atmospherically relevant conditions, only the lower volatility component is expected to partition appreciably to the aerosol phase.
"""

#dibutyl phtalate partition coefficient

	Mdbp = 278.348 #unit = g/mol" #Molecular weight
	ΔᵥₐₚH = 91.7 #kJ/mol
	p_dbp = exp.(ΔᵥₐₚH * 1000 / R .* (1 / 298.15 .- 1 ./ T)) #vapor pressure in P
	
	Kₒₘ_dbp = zeros(18)
	for w in 1:length(T), e in 1:length(p_dbp)
		Kₒₘ_dbp[w] = R * T[w] / (Mdbp * 1 * p_dbp[e] * 1000000)
	end
	Kₒₘ_dbp #m3/ug
	
#nonadecane

	#nonadecane

	Mnd = 268.5 #unit = g/mol" #Molecular weight
	ΔᵥₐₚH₂ = 95.8 #kJ/mol
	p_nd = zeros(18)
		p_nd = exp.(ΔᵥₐₚH₂*1000/R .* (1/298.15 .- 1 ./ T)) #vapor pressure in Pa
	p_nd
	
	Kₒₘ_nd = zeros(18)
	for n in 1:length(T), j in 1:length(p_nd)
		Kₒₘ_nd[n] = R * T[n] / (Mnd * 1 * p_nd[j] * 1000000)
	end
	Kₒₘ_nd #m3/ug
	
plot(T, Kₒₘ, c=:blue, lab = :none, xlabel="Kₒₘ", ylabel="Temperature")
	plot!(T, Kₒₘ_dbp, c=:red, lab = :none, xlabel="Kₒₘ", ylabel="Temperature")
	plot!(T, Kₒₘ_nd, c=:orange, lab = :none, xlabel="Kₒₘ", ylabel="Temperature")
	
md """
As we see by the plot the trends of Kₒₘ is pretty similar for all three species, with the bigger values for semivolatile organic compounds, which is consistent with the results of the paper. For the oxidation forms of the intermidiate volatile organic compounds (naphtalene used as a IVOC in the paper) increases significantly over the 10^9 times.

If we try to reproduce the same process with the saturation concentration C', which recalls the equation 2 from Pye et al. (2010), we will get following plot:
"""

Kₒₘₛₐₜ = 1 ./C_sat

plot(species,Kₒₘₛₐₜ, seriestype = :scatter, c=:blue, lab = :none, xlabel="Kₒₘ")

#I included the suggestion to see the correlation of the values obtained for each combination of temperature and pressure, and created a heatmap, but I considered it makes more sense to keep for loop since the vapor pressure is different and unique for each temperature.

Kₒₘ2 = R .* Tp[:,1] ./ (Mᵢ[1] .* γᵢ[1] .* Tp[:,2]' * 1000)

heatmap(Kₒₘ2,
	xlabel="Temperature", ylabel="Vapor pressure", 
		colorbar_title="Equilibrium partition coefficient",
		title="Kₒₘ dependance on Temperature and vapor pressure")
		
#The calculations provided help to understand the processes occurring in ambient atmosphere, equilibrium partitioning coefficient was calculated for different species with the dependence on temperature. It can be further used for modeling partition of low volatile compounds based on the data of the emissions of the low volatile compounds. However, GEOSChem model is already existing and was used in the Pye et al.(2010), but the estimates of the emissions and the contribution of each species is still not accurate enough, so more precise modeling could be conducted.
		
end
